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Half a century ago a reaction-diffusion system of two chemicals was introduced by Alan Turing to account for 
morphogenesis, i.e., the development of patterns, shapes and structures found in nature. Here we will discuss 
the formation of patterns and structures obtained through numerical simulation of the Turing mechanism in two 
and three dimensions. The forming patterns are found to depend strongly on the initial and boundary conditions 
as well as system parameters, showing a rich variety of patterns, e.g. stripes and spots (2D), and lamellae and 
spherical droplets (3D) arranged in structures of high symmetry, with or without defects or distortions. 


1 Introduction 


In the quest for understanding biological growth, it was 
Alan Turing (one of the key scientists of 20th century) who 
first demonstrated how a simple model system of coupled 
reaction-diffusion equations could give rise to spatial pat- 
terns in chemical concentrations through a process of chem- 
ical instability [1]. He showed that this kind of system 
may have a homogeneous stationary state which is unsta- 
ble against perturbations, such that any random deviation 
from the stationary state leads through diffusion to a sym- 
metry break. This process is called diffusion-driven insta- 
bility. Since complex spatial patterns are commonly found 
in nature [2], for example in animal skins and also in some 
polymer systems [3], it is quite natural to think that such pat- 
tern formations could be caused by some general physico- 
chemical process. Inspired by this notion, Turing systems 
have been proposed to account for pattern formation in var- 
ious biological systems [4], e.g. patterns on fish [5, 6], but- 
terflies [7] and lady beetles [8]. Furthermore, the first exper- 
imental evidence of a Turing structure was observed quite 
recently in a single-phase open chemical reactor [9]. 

The forms and variability of patterns generated by Tur- 
ing systems have been studied by assuming inhomogeneous 
diffusion coefficients [10], or by introducing domain curva- 
ture [11] or growth [12]. In addition, the Turing systems 
have been investigated [10, 13] from the point of view of 
symmetry breaking due to its potential relevance in under- 
standing the mechanisms responsible for symmetries found 
in nature. Sofar most of the computational modelling work 
on Turing systems has been done in two-dimensions (2D), 
while three-dimensional (3D) systems have remained little 
studied, as pointed out by De Wit et al. [14]. They investi- 
gated the so called Brusselator reaction-diffusion model in 
3D and demonstrated the possibility of additional high sym- 
metry structures, i.e. body centered cubic (BCC) and hexag- 
onally packed cylinders (HPC), and also a distorted lamel- 
lar structure which appears as a twist grain boundary. This 


latter structure embeds a minimal Scherk surface. Thus it 
is evident that 3D systems show a much richer set of mor- 
phologies than 2D systems. 


Recently we have also studied the effect of dimensional- 
ity by simulating a three-dimensional Turing system, which 
clearly shows an increased complexity in pattern formation 
[15]. While in 2D we obtain structures with spots organized 
hexagonally or with stripes, or some type of labyrinthine 
patterns, in 3D complex lamellar, spherical droplet patterns 
or their combinations are seen. We have also observed that 
initial conditions play a crucial role as to which stationary 
or long living pattern the system evolves. For example with 
certain initial introduction of one of the morphogens the sys- 
tem seems to evolve to a hexagonally arranged spherical 
droplet pattern or to a perfect lamellar structure, depend- 
ing on the parameters of the model system. Also we have 
studied situations in which the choice of system parameters 
indicate competition between spherical droplet and lamel- 
lar structures, and they give rise to complex hybrid patterns. 
In addition we have studied connectivity and clustering in 
the patterns [16] and the effect of noise on Turing structures 
[17]. 


These findings are in line with the observations by De 
Wit et al. [14] that 3D reaction-diffusion systems can show a 
wide variety of possible perfect and distorted patterns. How- 
ever, much remains to be explored and studied in more de- 
tail. Thus the purpose of this study is to explore the mor- 
phologies of patterns within the framework of generic Tur- 
ing system, in which it is possible to control pattern forma- 
tion with the parameters of the system. Also the effects of 
initial and boundary conditions are going to be investigated. 


In the next section we briefly present the generic Turing 
model and discuss its characteristics from the point of view 
of mode instability. Then we present the results of numeri- 
cal simulations in both 2D and 3D systems. Finally we draw 
conclusions. 
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2 Generic Turing Model 


As mentioned above a Turing system describes temporal be- 
haviour of the concentrations of two reacting and diffusing 
chemicals or morphogens. This can be represented in gen- 
eral by the following coupled reaction-diffusion equations 


U = 
V = 


DyV?U + f(U,V) 
DyV?V +g(U,V), (1) 


where U = U(#,t) and V = V(#,t) are the morphogen 
concentrations, and Dy and Dy the corresponding diffusion 
coefficients setting the time scales for diffusion. The reac- 
tion kinetics is described by the two non-linear functions f 
and g. 

In this paper we use the generic Turing model introduced 
by Barrio et al. [6], in which the reaction kinetics was devel- 
oped by Taylor expanding the non-linear functions around 
a stationary solution (Uc, Ve), defined by f(U.,V.) = 
g(Uc, Ve) = 0. If terms of the fourth order and higher are 
neglected, the reaction-diffusion equations can be written as 


D65V7u + au(l — riv?) + v(1 — rgu) 
v = 6V?vtv(B+ariuv)+uly+rev), (2) 


Ut = 


where u = U — U. and v = V — Y.. The parameters rı and 
T2 set the amplitudes of the non-linear cubic and quadratic 
terms, respectively, serving as control for favouring a certain 
pattern formation. In fact by gradually changing the non- 
linear parameters we observe a transition from spotty (2D) 
or spherical droplet (3D) patterns to striped (2D) or lamellar 
(3D) patterns. In the equation the quantity D is the ratio of 
the diffusion coefficients of the two chemicals, and 6 acts as 
a scaling factor fixing the size of the system. In the analysis 
of the model system we can without loss of generality fix 
a = —y in order to have only one stationary state in the sys- 
tem. For details about the modes of instability and the linear 
stability analysis of the model we refer the reader to Barrio 
et al. [6, 15]. Although the conditions for diffusion-driven 
instability to occur are widely known [4], it should be men- 
tioned that setting D Æ 1 is a necessary but not sufficient 
condition for the diffusion-driven instability in 2D and 3D 
systems [18], and that by restricting the parameter selection 
such that a € (0,1) and 6 € (—1,0) one is left with only 
two additional conditions, namely 


-1<B<-a, 


a — 2y Dj|a| > BD. 
Now from the linear stability analysis [6, 15] we obtain 
the following dispersion relation 


A’ ~(a+8—6k? (1+ D))A+(a—dDk?)(B—5k?) +a = 0. 
(3) 
Based on this equation one can find a region in k-space with 
positive growth rate, i.e., the eigenmodes u = uge and 
v = voe™ with eigenvalues A(k) > 0. In addition, we can 
analytically derive the modulus of the critical wave vector 


For the numerical integration of Eq. (2) the (finite) model 
system needs to be discretized. Then in a three-dimensional 
cubic system, the wave number is of the form 


= 2 
JB = and tn +n, © 


where L is the system size and nz, ny, nz are the wave num- 
ber indices (Note that in a two-dimensional system n, = 0). 
By adjusting the parameters and allowing only a few unsta- 
ble modes, one can end up with several different parame- 
ter sets. As in our earlier work [15] we chose the param- 
eters D = 0.516, a = —y = 0.899, 86 = —0.91 and 
ô = 2 corresponding to a critical wave vector ke = 0.45, 
and D = 0.122, a = 0.398, 8 = —0.4 and ô = 2 corre- 
sponding to k, = 0.84. 


3 Simulation Results 


For numerical simulations the spatial dimensions of the sys- 
tem were first discretized into a square or cubic lattice mesh 
with lattice constants dx = dy = dz = 1.0. The equa- 
tions of motion, Eq. (2) were iterated in time using the Euler 
scheme with time step dt = 0.05. Both zero-flux and peri- 
odic boundary conditions were used. 


In Fig. 1 we present the results for the 2D Turing sys- 
tem. In the left column we show the simulation results of 
Turing patterns for four sets of non-linear parameters (71 
and rz in Eq. (2)). By slowly changing the strength of the 
quadratic nonlinear interaction term (r2) one can observe 
transition from stripes to spots. Our preliminary results sug- 
gest that in the limit of infinite simulation time the transi- 
tion between stripes and spots could become discontinuous. 
However, with a finite simulation time one always observes 
intermediate transient states (Figs. 3B and C) for fixed crit- 
ical amplitudes of the nonlinearities. In the proximity of 
these critical values one observes critical slowing down of 
the dynamics: the closer one is to the critical values of the 
nonlinear terms the longer it takes for the system to reach 
either regular striped or spotty morphology. 


One can get further insight into the transition by taking 
the discrete Fourier transform of the concentration data into 
the wave vector space, k, i.e., 


ak) => oe. 6) 


This approach has been used earlier, e.g. for reaction- 
diffusion systems [19], Turing patterns [20], and to char- 
acterize the evolution of patterns [21]. The quantity |(k)| 
corresponds to a diffraction pattern. The results are shown 
in the right column of Fig. 1, which shows a sequence of 
diffraction patterns corresponding to the concentration space 
Turing patterns on the left. 
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Figure 1. 2D Turing system: Transition from stripes to spots in the 
real space (left column) and the corresponding diffraction patterns 
in wave vector (kz, ky)-space (right column). The patterns have 
been obtained from random initial conditions after 500000 itera- 
tions in a 200 x 200 system with ke = 0.45. The nonlinear param- 
eters were rı = 2.0 and A) r2 = 0, B) r2 = 0.162, C) r2 = 0.167 
and D) re = 0.367. In the real space black corresponds to areas 
dominated by chemical U and the lighter color by chemical V. In 
the reciprocal space (kz, ky) = (0,0) is in the middle. 


In Fig. 1A, the diffraction intensity is to k,,-direction be- 
cause the stripes are aligned in y-direction. The distance of 
these diffraction peaks from the origin gives the length of 
the wave vector of the unstable mode, while its width in the 
perpendicular direction describes deviations of stripes from 
the principal directions due to stripe being tilted and bent. 
In Fig. 1B the diffraction intensity becomes more spread 
around the diagonal in (k,k,)-space, as a result of nucle- 
ating defects and spots. After that in Fig. 1C the diffrac- 
tion intensity starts splitting into separate peaks due to more 
spots forming, then developing into six separate equidistant 
(from the origin and from each other) intensity peaks as ev- 
ident in Fig. 1D. This hexagonal symmetry in the reciprocal 
space is because the system evolves towards regular spotty 
pattern with predominantly triangular symmetry in the real 
concentration space. However, these diffraction peaks are 
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Figure 2. 3D Turing structures: Spherical droplets arranged in 
BCC symmetry (ke = 0.45, rı = 0.02, r2 = 0.2), planar lamel- 
lar structure (ke = 0.84, rı = 3.5, ro = 0) and HPC structure 
(ke = 0.84, rı = 0.25, r2 = 0.15). 


somewhat spread, which indicates that the triangular sym- 
metry is not perfect over the whole system. 

In three dimensions the pattern selection is not as 
straightforward as in two dimensions. Instead of perfectly 
organized spherical structures and planar lamellae, one may 
obtain structures that establish order in more general ways. 
In the case of spherical droplet structures, the droplets are 
not always organized on a perfect lattice. The lamellae are 
rarely planar, but more complex minimal surfaces or twisted 
grain boundaries arise, c.f. [14]. In addition, the intermedi- 
ate structures are not separate domains of spherical droplet 
structures and lamellae as spots and stripes in 2D, but a va- 
riety of stable structures is obtained. For typical 3D Turing 
structures we refer the reader to [15]. 

Some of the 3D structures generated by the generic Tur- 
ing model resemble those that have been predicted and simu- 
lated in diblock copolymer melts [22, 23] and also those ob- 
tained by De Wit et al. [14] using the Brusselator model. In 
Fig. 2 we show specifically three Turing structures that are 
also observed in polymer melt systems: spherical droplets 
organized into BCC lattice, planar lamellar structure, and 
HPC structure. These structures were obtained by numeri- 
cal simulation of the generic Turing model of Eq. (2) 
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Figure 3. Minimal isoconcentration surface for connecting two 
crosswise planar structures. The close-up reveals the morphology 
of the Scherk surface. 


after 400000 iterations in a 50 x 50 x 50 system. The ini- 
tial conditions and the system parameters had to be adjusted 
carefully to obtain these structures. However, the structures 
proved to be stable against noise once established. 

For the BCC structure in Fig. 2 the initial concentration 
field had the BCC symmetry, which made the system find 
the BCC structure of larger spherical droplets. This structure 
we could not, however, get by starting from a random initial 
configuration of morphogens. For the planar lamellar struc- 
ture, chemical U was introduced only in the midplane of 
the system and V was distributed uniformly over the whole 
system. For the same set of parameters but starting from 
random introduction of both chemicals we find a lamellar 
minimal surface, which is characterised to have both posi- 
tive and negative local principal curvature. As for the HPC 
structure, it was obtained starting the simulation from ran- 
dom introduction of both chemicals. The used parameter 
values can be found from the caption of Fig. 2. 

There is an optimal minimal surface for connecting per- 
pendicular planes such that the local principal curvatures (c1 
and c2) have opposite signs yielding zero local mean curva- 
ture (H = (cı + c2)/2 = 0). This surface is called Scherk 
surface, which is a well-known minimal surface solution for 
twisted lamellar surfaces. The stability of a Scherk surface 
has previously been observed in the case of the Brussela- 
tor model [14] by initializing the minimal surface into the 
system. In Fig. 3 we show the final structure obtained as 
a result of a simulation run of a system initializing one of 
the chemicals into two planar domains set perpendicularly. 
Many twisted connections of lamellar surfaces are found. 


Based on the close-up of one of the connections we can ob- 
serve that also spontaneous formation of the Scherk surface 
is possible in the generic Turing system. The parameters in 
this case were rı = 3.5 and r2 = 0 for mode ke = 0.84. 


4 Summary 


In summary we have studied the pattern formation in a 
generic Turing model in two and three dimensions and under 
different initial and boundary condition. The generic Turing 
model is computationally very interesting and rich because 
one can change the nonlinear parameters to control the pat- 
tern formation and explore the space of pattern the system 
produces. In 2D we have observed striped and spotty pat- 
terns and a transition between them when the nonlinear pa- 
rameters are changed. This transition turns out to be very 
sharp, or possibly even discontinuous at the limit of infi- 
nite time. In 3D the number of possible stationary patterns 
turned out to be much larger than in 2D. We have observed 
spherical droplet pattern (organized either somewhat irregu- 
larly or in hexagonal arrangement), regular BCC pattern of 
spherical droplets, regular HPC structure, perfect lamellar 
structure, distorted lamellar structures which carry charac- 
teristics of minimal surface, and Scherk surface structure of 
minimal surface solution for twisted lamellar surfaces. 

To conclude it is evident that Turing systems have very 
general ways of establishing order since we have observed 
spontaneous formation of a wide variety of stable minimal 
surfaces in addition to the Scherk surface. All these surfaces 
have the characteristic length fixed by the parameters and 
predicted by linear analysis, but the morphologies can be 
very different. This interesting problem will be addressed in 
more detail in future research. 
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